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The Spalart-Allmaras and the Menter SST k- u turbulence models are shown to have the undesirable char- 
acteristic that, for fully turbulent computations, a transition region can occur whose extent varies with grid 
density. Extremely fine two-dimensional grids over the front portion of an airfoil are used to demonstrate 
the effect. As the grid density is increased, the laminar region near the nose becomes larger. In the Spalart- 
Allmaras model this behavior is due to convergence to a laminar-behavior fixed point that occurs in practice 
when freestream turbulence is below some threshold. It is the result of a feature purposefully added to the 
original model in conjunction with a special trip function. This degenerate fixed point can also cause non- 
uniqueness regarding where transition initiates on a given grid. Consistent fully turbulent results can easily 
be achieved by either using a higher freestream turbulence level or by making a simple change to one of the 
model constants. Two-equation k-u> models, including the SST model, exhibit strong sensitivity to numerical 
resolution near the area where turbulence initiates. Thus, inconsistent apparent transition behavior with grid 
refinement in this case does not appear to stem from the presence of a degenerate fixed point. Rather, it is a 
fundamental property of the k-uj model itself, and is not easily remedied. 


I. Introduction 

Within the aerodynamics community, the Spalart-Allmaras (SA) 1 and the Menter shear stress transport (SST) 2 k-us 
turbulence models have been widely-used and trusted for Reynolds-averaged Navier-Stokes (RANS) computations for 
over a decade. The SA model is a one-equation model based primarily on empiricism and on dimensional analysis 
arguments. It works well for both attached and for many separated flows. It was one of the first field-equation models 
to gain wide acceptance into modem aerodynamic application CFD codes. The SST model combined the robustness 
of the k-uj turbulence model 3 near walls with the capabilities of the k-e model 4 away from walls (also removing 
the sensitivity to variations in farfield lj boundaiy conditions 5 ). It also included the effects of turbulent shear stress 
transport through Bradshaw’s assumption 6 that the shear stress in a boundary layer is proportional to the turbulent 
kinetic energy. This “SST” modification was the key that gave the model its ability to predict separated flows much 
better than the standard k-uj form. Both the SA and SST models have been remarkably successful and stable (very few 
changes proposed or required) over the years. 

Recently, Rumsey et al. 7 showed that certain forms of the low Reynolds number two-equation k-e model could 
exhibit apparent transition behavior that is dependent on initial conditions and method of solution (such as different 
procedures of mesh sequencing). They determined the numerical reasons behind the long-recognized inconsistent 
transition prediction behavior of the k-e model (see, e.g., Abid 8 ). The apparent laminar flow region predicted by the 
model upstream of transition was termed “pseudo-laminar” because, although the eddy viscosity was predicted to be 
near-zero (yielding laminar results for all intents and purposes), the behavior of the near-zero turbulent quantities were 
incorrect relative to each other in the laminar limit. The pseudo-laminar state could be reached in some regions of the 
flow solely because k = e = 0 was a stable solution to the equations under certain conditions, and not because of any 
physics built into the model. 

CFD practitioners are generally aware that most transport turbulence models “transition” on their own. 9-11 The 
location of this apparent transition generally has no relationship with experimental results. In other words, in computed 
results there is usually a region of laminar flow followed by turbulence, even if a CFD code is run in fully-turbulent 
mode, with the turbulence equations active everywhere in the flow field. For the SA and SST models at reasonably 
high Reynolds numbers, this apparent transition location tends to occur very far upstream in most applications, and 
the models seem to avoid the aforementioned problems often seen with the k-e model. 
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However, in light of the recent analysis of the k-e model, 7 it is appropriate to apply a similar analysis to the SA 
and SST models in an effort to shed some light on their ability to produce undesirable apparent transition behavior. 
The main question to try to answer is whether there are circumstances for which the SA and SST models can yield 
inconsistent results in terms of the extent of laminar flow predicted. If present, such an inconsistency can have negative 
consequences when attempting to draw conclusions from computed results. 

II. Analysis of the Turbulence Models 


A. Spalart-AIlmaras 

The one-equation SA model is written in terms of the turbulence quantity 0. 
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The eddy viscosity is given by v t — 0f v \ , ft is the magnitude of the vorticity, d is the distance to the nearest wall, and 
the constants are: Cm = 0.1355, C b2 = 0.622, h = 0.41, a = 2/3, C ts = 1.2, C t4 = 0.5, C w2 = 0.3, C u>3 = 2, 
C vl =7.1, and C w \ — C b \/t\~ + (1 -b O b 2 )/cr. 

As noted in Spalart and Allmaras, 1 the f t2 term was an ad hoc numerical fix that makes 0 = 0 a stable solution 
to the equations with a small basin of attraction (also note that its constants C t2 and Cu changed values between the 
AIAA paper reference and the subsequent journal article reference). The analysis below confirms this characteristic. 
The reason the model developers originally wanted this behavior was due to their use of a “trip function,” /, x At/, not 
described above but given in the original reference. Because the SA model without f r2 would often trip to turbulence 
upstream of where the numerical trip was set for transitional flows, Spalart and Allmaras invented the f r2 term to 
delay it. However, arguably most present-day users of SA do not employ the trip function, but rather run the model in 
fully-turbulent mode (or else force desired laminar regions by zeroing out the production term). In the Results section, 
potential undesirable consequences of the presence of the f r > term for fully-turbulent computations will be shown. 

For the first part of the analysis here, the homogeneous form of the one-equation model is used. As discussed in 
Rumsey et al., 7 this is a useful approximation because the log-layer or equilibrium layer of a turbulent boundary layer 
has the same dynamical characteristics as a homogeneous flow. The dynamic similarity is exploited by treating the 
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flow as “locally homogeneous.” In other words, although the mean shear varies throughout a boundary layer, in the 
analysis the assumption is made of a locally fixed mean shear at a given location. 

The homogeneous form of the SA equation is 


where 
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By inspection, it is clearly seen that 0 = 0 is a critical point that satisfies the equation in the steady-state. The stability 
of this solution is investigated by solving for dR/dv in the vicinity of v = 0: 
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Because C t3 = 1.2, Eq. (13) is always negative, which means that the null critical point is always stable. Therefore, 
whenever v gets in the vicinity of v = 0, the solution to the partial differential equation may be attracted toward v = 0. 

This stability can be demonstrated numerically by solving for the instantaneous time rate-of-change of the tur- 
bulence quantity for a wide array of possible current state values. This is done by first writing Eqs. (10) - (11) in 
nondimensional form: 
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(The quantity r also requires a (M/Re) factor: r = (M / Re)i>' / (S' k 2 d! 2 ) .) In these expressions, v f = z>/V r ef> 
t! = £a re f/L re f, d! = d/L re f , and Q' ~ ^L re f/a re f (with a re r representing the refe rence speed of sound). 

Next, a variety of values are input for the independent variables d f yjRe/M, and iY . An example solution map 
is shown in Fig. 1. Here, setting Q f = 1000, a contour plot of dv* jdt! is made for an array of v 9 and d' yjRe/M 
values. The solution map shows that dv 1 jdt 1 < 0 for the region below and to the left of the neutral line (the solid 
line in the figure). This means that for an initial condition in this region the solution variable v’ is driven lower as the 
numerical iteration procedure progresses. Above and to the right of the neutral line, dv 1 jdt! > 0, which means that v’ 
is driven higher. Thus, for d! yjRe/M to the right of the dot in the figure, there are two possible stable solutions. The 
portion of the neutral line above the dot represents the turbulent solution to the homogeneous form of the SA equation. 
For a given d! \J Re/ M > 0.09 or so, the solution to the equation will converge to the corresponding />' on the stable 
neutral line, provided that the initial condition is somewhere above the unstable neutral line. Otherwise, the solution 
will converge to the degenerate laminar v 1 = 0 v alue. 

Although not shown, for larger d 9 y/Re/M, the lower part of the neutral line approaches a constant value near 
v’ » 0.6, and the upper part of the neutral line continues upward and to the right. Also, for different values of the 
general character of the solution map remains the same as that shown in the figure, but the specific contour shapes and 
levels shift. 

This exercise using the homogeneous form of the SA equation is instructive, but not realistic. In the full form of 
the equation, advection and diffusion are locally active and play a significant role. These terms can have basically three 
effects on the solution map, when added to the right-hand-side of Eq. (14). First, if the sum of the terms is negative. 
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then the general character of the solution map remains unchanged from the homogeneous map: only its levels and 
location of the neutral lines are different. Second, if the sum of the terms is positive and relatively large compared to 
the other terms in the equation, then the possible problem of laminar solution behavior goes away. An example solution 
map is shown in Fig. 2. For this map. ft' is again taken as 1000, but this time an assumed nondimensional advection 
plus diffusion contribution is taken to be a constant value of 10, for demonstration purposes. The map shows that in 
this case there is only one neutral line; it is stable and corresponds to the turbulent solution. Any initial condition will 
converge to this turbulent solution. For fully turbulent flow, this situation is the desired one: the numerics are such that 
the solution always converges to the intended, unique value. 

Third, very interesting non-unique behavior occurs if the sum of the right -hand-side advection plus diffusion is 
positive but not too large compared to the other terms. An example solution map is shown in Fig. 3. For demonstration 
purposes, ft' is again taken as 1000, and nondimensional advection plus diffusion is taken to be a constant value of 3. 
In this case, there are three distinct neutral lines. The lower line oriented horizontally is stable, and will be reached for 
any initial condition below the unstable neutral line. This lower stable neutral line has very low levels of i>\ yielding 
even smaller levels of eddy viscosity. Thus, convergence to this line represents a solution that behaves like laminar 
flow. It will be shown in the Results section that this solution behavior can cause the apparent transition location for a 
fully turbulent computation to vary noticeably with simple grid refinement. 

B. Menter SST k- ^ 

The two-equation SST model is written as follows: 
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and 7, = 0.55317, 72 = 0.44035, a k i = 1.17647, a k2 = 1, cr^i = 2, a u2 = 1.16822, ft = 0.075, ft = 0.0828, 

k = 0.41 , ai = 0.31 , 0* = 0.09, fi is the magnitude of the vorticity, and d is the distance to the nearest wall. Because 

many CFD codes use the approximation that P k = yu f £1 2 (as recommended by Menter 12 ), it is employed here as well. 

As an initial step, we look briefly at the homogeneous form of the SST model. Because we are concerned primarily 
with transition and the inner part of the boundary layer, we perform this initial analysis with F\ = F 2 = 1. Therefore, 
the cross-diffusion term in Eq. (17) is zero. Nondimensionalizing the turbulence quantities via: 

k' = k/a 2 . ef (29) 

J = ujn ie {/{p ret a%f) ( 30 ) 

and then, for convenience, further normalizing these variables, along with time, via: 
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the homogeneous form of SST inside the boundary layer can be written: 
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Similar to the SA model, a degenerate solution satisfies these equations in the steady state: k * — 0 9 UJ* — yj 7i / . 

This value of u;* is less than \/a\ for the values of the given constants, so a\k* < k* /w*, and the stability of the 

degenerate solution can be found from the eigenvalues of the coefficient matrix of the linear system (Eqs. (34) and 
(35)): 
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(36) 


Plugging in the degenerate solution, the eigenvalues are found to be real with one positive and one negative: this 
indicates that the solution is a saddle point. 7 Generally, a saddle point is not stable in practice because orbits approach 
the critical point along one eigenvector, but then recede along the eigenvector associated with the unstable solution. 
Thus, this analysis suggests that the SST model should not be attracted to the degenerate solution. 

However, as with the SA model, it is important to also consider the effects of the advective and diffusive terms in 
the SST model in order to truly begin to understand its dynamical behavior. The easiest way to do this is to perform an 
actual computation, and monitor the values of all the terms as the solution progresses. This has been done, and detailed 
results will be shown in the Results section. However, for the purposes of analysis and attempting to get a picture of the 
overall k*-u> * solution space, we make the assumption that the dominant turbulent transport effects act over a distance 
that represents the length scale characterizing the large turbulent eddies, 13 given by i = k s ^ 2 /e = k 1 ^ 2 /uj. Then, the 
nondimensional advection and diffusion terms in the turbulent region can be approximated by: 
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u>— eqn cross diffusion : 2(1 - F X )~ P . ^L^L ~ O ( ~ O (w 12 ) = (C^.cd)^' 2 ( 41 ) 

dx'j jdij \° J 

where fi\ = Pt//Xref> and C k . A , CL.d, CL v 4, CL,d, and CL.cd are (unknown) coefficients that can be determined 
at instantaneous locations and times in a particular flow field via numerical computations. Also, Uj represents the 
(nondimensional) mean flow velocity components, which are taken as locally constant and are absorbed into the 
coefficients C k . A and CL,a- Employing the normalizations from Eqs. (31) - (33), the inhomogeneous equations can 
be approximated as: 
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(with F\ and F 2 no longer assumed to be 1). These equations are admittedly very crude, because of the approximations 

made in Eqs, (37) - (41). In fact, analysis of the actual advection and diffusion terms from CFD results at several 

locations in the bounday layer confirmed that the coefficients C k , A , etc., do not come out to be constant. However, 
the length scale / = k : ^ 2 /s yielded a fairly reasonable approximation that was much better than other choices. The 
solution to Eqs. (42) and (43) can be found. It is the intersection of the A~*-nullcline, given by: 
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and the u*-nullcline, given by: 


\[3-(C^,a) k* 1/2 — CL,d - C^.cd) 

The analytic form of the eigenvalues associated with this solution is very complicated, and does not yield a clear 
conclusion regarding the stability of the solution. However, it is easy to numerically find the eigenvalues for a point 
in the turbulent boundaiy layer of a typical computation. In all cases investigated to date, the solution has turned 
out to be a saddle point. An example phase-plane portrait is shown in Fig. 4. This figure shows the &*- and u;*- 
nullclines using coefficients obtained at a point in a converged turbulent boundary layer where nondimensional eddy 
viscosity fj' t = 26.47524. In this flow field, Re = 1 x 1() 7 and M = 0.3. The converged solution at this point 
(circled in the figure) has a value of k f — 0.7312424 x 1()~ 3 , u/ = 0.2386361 x 10~ 4 , and nondimensional vorticity 
magnitude O' = 265.3249. Thus k* = 0.1940168 and uj* = 2.998036. Other coefficients obtained from the numerical 
solution were: C kJ) = -0.3368563 x 1(T 2 , C ktA = -0.4419278 x 10“ 2 , CL,p + C^cd = 0.9705072 x Hr 2 , 
C^a = 0.1760321 x 10 -2 , F\ = 0.9879624, and — 0.9999944. The figure also shows the local trajectories 
(defined by dk*/dt* and c)u;*/dt*) of the turbulence variables in phase-space. 

It is clear from the figure that the solution is a saddle point: solution trajectories approach from above and below, 
but recede in both directions along the u;*-nullcline. This is somewhat surprising, because at a saddle point an orbit 
is usually attracted at first but repelled later on. This would seem to imply that the SST model is not strictly stable. 
However, nonlinear and nonlocal effects were not accounted for in the analysis, and it turns out that these act to 
maintain stability at the solution point. For example, using a numerical experiment from the example in Fig. 4, when 
k* at the solution point was perturbed slightly to the left, then the local spatial k* derivatives also changed, yielding 
different local advective and diffusive terms which temporarily shifted the intersection point of the nullclines even 
further to the left; as a result the solution was driven back to the right, toward its original position. In any case, SST 
certainly is numerically stable in practice. 

Although not shown here, the standard A -u/ model 3 yields a phase-plane portrait very similar to that of SST, with 
its solution also appearing as a saddle point. Thus, its behavior is expected to be similar. 

Now the following question arises: what are the effects of this model’s seemingly tenuous linearized-solution 
stability on its apparent transition behavior? The answer is not entirely clear from this analysis, but in the Results 
section some details will be shown that will hopefully shed some light on this question. 
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III. Numerical Method 


The computer code CFL3D 14 solves the three-dimensional, time-dependent compressible RANS equations with 
an upwind finite-volume formulation (it can also be exercised in two-dimensional mode of operation for 2-D cases). 
Upwind-biased third-order spatial differencing is used for the inviscid terms, and viscous terms are centrally differ- 
enced. The code originally solved the thin-layer form of the equations (in each coordinate direction), but the full 
Navier-Stokes terms (i.e., cross-derivative terms) have recently been added. All solutions shown below use the full 
Navier-Stokes terms, although thin-layer has also been employed and makes almost no difference for the particular 
cases run. 

The CFL3D code is advanced in time with an implicit approximate factorization method. The implicit derivatives 
are written as spatially first-order accurate, which results in block tridiagonal inversions for each sweep. However, 
for solutions that utilize Roe flux-difference splitting, 15 the block tridiagonal inversions are further simplified using a 
diagonal algorithm with a spectral radius scaling of the viscous terms. 

The turbulence models, including SA and SST, are solved uncoupled from the mean flow equations using im- 
plicit approximate factorization. Their advective terms can be solved using either first-order or second-order upwind 
differencing, with first-order the default for the code. 


IV. Results 

In order to demonstrate inconsistent apparent transition behavior of the SA and SST models, it was necessary 
to perform computations on very fine grids. Easily noticeable inconsistencies did not show up using coarser grids. 
The model problem chosen was the computation over the front upper 45% of a NACA 0012 airfoil, at M = 0.3, 
Re = 10 7 based on total chord length, and a = 0°. Fully-turbulent computations were performed, which means 
that the turbulence model equations were active everywhere in the flow field. The freestream eddy viscosity was set 
using a typical value 16 of turbulence Reynolds number Hr, oc = kto/i^o ©£oo) — 0.1, yielding /^ <00 = 0.009. For the 
SA model, this was achieved using v ^ = 1.341946. It will be shown below that for the SA model, the freestream 
turbulence level chosen can be an important consideration. For the SST model, the freestream turbulent quantities 
used were: k f OQ = 9 x 10 -9 and uo'^ = 1 x 1CT 6 . 

The finest grid was 769 x 1025, with coarser grids created by successively removing every other grid point from 
the next finest level. Thus the coarser grid sizes were: 385 x 513, 193 x 257, 97 x 129, and 49 x 65. The grids had 
a farfield extent of 50 chords, and the minimum normal spacing at the wall yielded an average y - b level of 0.05 for 
the 769 x 1025 grid, 0.1 for the 385 x 513 grid, 0.2 for the 193 x 257 grid, 0.4 for the 97 x 129 grid, and 0.8 for the 
49 x 65 grid. 

A picture of the coarsest grid is shown in Fig. 5. Note that this coarsest level is on a par with typical grid resolution 
used for many turbulent CFD computations over airfoils and wings. For boundary conditions, the part of the grid 
coming into the nose of the airfoil was set as a symmetry plane, the airfoil body used no-slip adiabatic wall conditions, 
the downstream exit plane used extrapolation, and the farfield (upstream and above the airfoil) used a farfield Riemann 
invariant condition. For most of the computations, each grid level was run by restarting from the previous coarse grid 
solution, as is typically done for full-approximation scheme (FAS) multigrid simulations. 17 All results were converged 
5-6 orders of magnitude on each grid (the final density residual dropped to 10 _1 ° - 10 -11 ). 

A. Spalart-Allmaras 

Fig. 6 shows computed surface skin friction coefficients on the airfoil using SA for each of the grids. As the grid was 
refined, each successive solution produced a larger “laminar region” prior to going turbulent. Only on the finest 769 x 
1025 grid did the transition prediction appear to stop changing significantly. For the coarsest grid, the nondimensional 
eddy viscosity first exceeded 1.0 in the boundary layer near xjc — 0.0044. For successively finer grids, the location 
was 0.0056, 0.0099, 0.0111, and 0.0112, respectively. As will be shown next, this apparent later transition on the 
finer grids is a numerical artifact related to the existence of the laminar-behavior stable neutral line for this turbulence 
model. 

Solution maps can be created for this realistic computation by monitoring the various independent quantities during 
the course of the solution’s progress. In this case, the quantities we re recor ded during the solution on the 193 x 257 
grid, at a point in the boundary layer near x/c = 0.00694, where d'y/RejM = 0.1157. At this position, the converged 
solution on the 97 x 129 grid was turbulent (i.e., the peak nondimensional eddy viscosity level - located somewhat 
further from the wall at this location - exceeded 1), but during iterations on the 193 x 257 grid this location became 
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laminar. 

Fig. 7 shows the descent of the v ! quantity at this location from near 3.9 initial value to a converged value of 0.205. 
(These correspond to nondimensional eddy viscosity levels of about 0.48 and 4 x 10~ 6 , respectively). This solution on 
the 193 x 257 was restarted from the converged solution on the 97 x 129 grid at multigrid cycle number 8001 . Next, by 
using the actual instantaneous values for O' and nondimensional advection plus diffusion at this location, the solution 
maps were plotted at multigrid cycle numbers 8940, 9120, and 11000 (at cycle 8940: v f — 3.1064, O' = 5365.3, and 
advection plus diffusion — -748.8; at cycle 9120: i>' = 2.313,0' = 5399.4, and advection plus diffusion = -1137.4; 
at cycle 1 1000: v* = 0.205, O' = 5500.3, and advection plus diffusion = 24.5). These maps are shown in Figs. 8, 9, 
and 10. At cycle 8940, the solution map shows that the location of the current value of V (indicated by the open circle) 
lies slightly above and to the left of the stable neutral line in the map space. Thus, its value is being driven lower with 
successive iterations. As its value decreases, the local advection and diffusion also change significantly such that the 
stable neutral line is driven further to the right (Fig. 9). Finally, at convergence (Fig. 10), the final value of />' lies on 
the stable neutral line representative of a laminar-behavior solution. 

Because two stable attractors exist, it is possible to obtain different solutions with SA depending on numerical 
factors such as initial condition (I.C.) and method of solution. An example is shown for the 97 x 129 grid in Fig. 11. 
The result labeled “I.C. #1” is the same result on this grid level shown earlier in Fig. 6. The result labeled “I.C. #2” 
was obtained from an initial condition using a converged result on the 193 x 257 grid interpolated onto the coarser 
level, then re-iterated to convergence. Two different apparent transition locations result: the first near x/c = 0.0056 
and the second near x/c = 0.0097. 

The numerically-induced laminar behavior of the SA model is easily avoided in two different ways. The first 
method, as suggested by Spalart, 10 is to employ higher levels of freestream turbulence. For example, using — 3 
(fi f toc = 0.2104), the problem of inconsistent transition was eliminated, as shown in Fig. 12. The apparent transition 
location remained near the airfoil leading edge (near x/c = 0.004 at this Reynolds number) for all grid levels. 

The earlier SA analysis also suggests a second method for avoiding the problem. By setting the constant C t 3 - 0, 
the right-hand side of Eq. ( 13) remains always positive, and the null critical point is unstable {y = 0 is not an attractor). 
The same NACA 0012 airfoil case was run with %(X> = 0.009 and C t:i = 0. Although not shown, results again show 
consistent behavior, and skin friction results appear identical to those shown in Fig. 12. 

B. Menter SST k-x 

The same NACA 0012 case was run using the SST model in fully turbulent mode. Results are shown in Fig. 13, and 
are similar to those using SA: as the grid was refined, the apparent transition location moved further aft. Even on 
the finest 769 x 1025 grid, the location still appeared to be influenced by the grid density. For the coarsest grid, the 
nondimensional eddy viscosity first exceeded 1.0 in the boundary layer near x/c — 0.0054. For successively finer 
grids, the location was 0.0062, 0.0076, 0.0086, and 0.0089, respectively. It should be noted that, unlike the SA model 
above and the k-e model reported in Rumsey et al., 7 SST did not appear to exhibit sensitivity to initial conditions or 
method of solution. In other words, a given grid always gave the same result. 

The dynamics of the SST solution behavior near the leading edge region was investigated by monitoring various 
quantities during the solution pr ogress o n the 193 x 257 grid, at a point in the boundary layer near x/c = 0.00694 
(distance from the wall was d'y/Rc/M = 0.2285). At this streamwise location, the solution on the 97 x 129 grid 
was turbulent (peak nondimensional eddy viscosity exceeded 1 in the boundary layer), but during iterations on the 
193 x 257 grid the peak eddy viscosity dropped below 1. 

Fig. 14 shows the changes in fi ' t , k * , and u/* through iteration to convergence. The value of a;* changed relatively 
little, but k* decreased by nearly 75%, causing a similar dramatic decrease in eddy viscosity. At multigrid cycle 
number 8500, the phase-plane portrait in Fig. 15 shows that the unconverged values of k* and ui* (indicated by the 
open circle) lie to the left of the intersection point of the nullclines. Thus, the solution in phase-space continues to be 
driven to the left. As it moves, the nullclines also change, but their intersection remains to the right of the solution until 
convergence is reached. Finally, at convergence (Fig. 16), the solution and the intersection of the nullclines coincide, 
as they must. 

Unlike SA, the SST model continued to exhibit inconsistent transition results even when the freestream turbulence 
level was raised to levels above 1. For example. Fig. 17 shows surface skin friction coefficients with fi' t ^ — 1.294 
( k ^ - 1.294 x 10 “). Compared with results in Fig. 13, these results using higher freestream turbulence show less 
grid influence, but the influence is still visible nonetheless. 

From these examples, it seems evident that the SST model does indeed exhibit a sensitivity of its predicted tran- 
sition location to grid density, but this sensitivity is not due to a null critical point acting as an additional attractor. 
Rather, there is only one possible solution, and the model dynamical behavior itself (near the area where turbulence 
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initiates) appears to be particularly sensitive to small changes induced by different numerical resolution. One of the 
reasons for this sensitivity may be the fact that, in a linear sense, the solution is a saddle point in k*-u>* phase-space. 
No simple cure was found to remedy this problem. It seems that it would be more desirable to have the system of 
turbulence equations exhibit strong stability, such as that shown in Fig. 18. Here, in this hypothetical plot (not based 
on any real model), the solution trajectories are all strongly attracted in toward the fixed point. 

As mentioned earlier, although the analysis and examples were done here specifically for SST, similar results are 
exhibited by the general form of the k-u; model as well. 


V. Conclusions 

In conclusion, both the Spalart-Allmaras and the Menter SST k- u turbulence models were shown to have the 
undesirable characteristic that, for fully turbulent computations, their apparent transition location can vary with grid 
density. As the grid was refined for an example airfoil computation, an apparent laminar region near the nose of the 
airfoil grew larger. In the SA model this behavior was due to a feature purposefully added to the original model in 
conjunction with a special trip function, which many people employing the model do not even use. Solution maps 
were used to show that this feature can cause convergence to a laminar-behaving fixed point. Because two stable fixed 
points exist, even a given grid can produce non-unique results regarding where transition initiates. It was also shown 
that consistent fully turbulent SA results can be achieved in practice by either using higher freestream turbulence levels 
(i/^ = 3 or greater), or by setting the constant C t 3 in the model to 0 . 

The SST model’s dynamical behavior (near the area where turbulence initiates) appears to be particularly sensitive 
to small changes induced by changes in numerical resolution. Thus, its apparent transition location change with grid 
refinement seems to be a fundamental property of the k-oj model itself, and is not easily remedied. Raising freestream 
turbulence levels did not fix the problem. Phase-plane portraits, which included the nullclines of the k-u; equations, 
were used to demonstrate the SST solution behavior. This model does not appear to suffer from the problem of possible 
non-uniqueness on a given grid. 

It is important to note that the problems illustrated in this paper do not appear to be very significant in practice. 
For example, the transition location movement appears to be fairly small for airfoils relative to the airfoil chord. In 
other words, generally these two turbulence models tend to transition near enough to the leading edge of airfoils and 
wings so that the overall solution remains globally consistent. However, this may be problem and/or Reynolds number 
dependent. In any case, running the SA model with > 3 for computations that do not use the trip function appears 
to be a good idea. Setting C t3 = 0 is also a viable choice, but may be less desirable from the point of view of model 
backward compatibility, because it is possible that the model behavior in and near separation could also be affected 
by this change. For k-u; models, including SST, users should be aware of the apparent transition-location dependence 
on grid density. In particular, as computers become more powerful and finer grids are more routinely employed, the 
problem behavior may manifest itself to a greater degree. 
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Figure 1. Homogeneous SA solution map of 80' /dt' contours for a range of initial conditions for w' and d! yj Re/M; the solid line 
represents dv' /dt 1 = 0; ft' = 1000. Initial values below and to the left of the neutral line are driven lower, and initial values above and to 
the right of the neutral line are driven higher. 
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Figure 2. Inhomogeneous SA solution map of dv' /dt' contours for a range of initial conditions for v' and d! yj Re/M ; the solid line 
represents dv' /dt' = 0; Q' = 1000; advection plus diffusion contribution taken to be 10. 
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Figure 3. Inhomogeneous SA solution map of dv' /dt' contours for a range of initial conditions for v* and d! yj Re/M ; the solid line 
represents dO' / dt' = 0; Q' = 1000; advection plus diffusion contribution taken to be 3. 
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Figure 4. Inhomogeneous phase-plane portrait of the SST model in k* and lj* space for a typical turbulent solution point (circled). 
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Figure 5. Portion of coarsest NACA 0012 grid, 49 X 65 . 
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Figure 6. Computed surface skin friction coefficient for the NACA 0012 airfoil on various grids, SA model with n' = 0.009. 



multigrid cycle number 


Figure 7. Convergence of v' at specific location in the boundary layer on 193 x 257 grid, SA model with n't oc = 0.009; solution maps will 
be plotted for the three times indicated in this figure. 
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Figure 8. SA solution map for NACA 0012 ease on 193 X 257 grid at a specific point in the boundary layer near x/c = 0.00694, at multigrid 
cycle number 8940. Open circle indicates the particular (unconverged) solution at this time. 



Figure 9. SA solution map for NACA 0012 case on 193 X 257 grid at a specific point in the boundary layer near x/c = 0.00694, at multigrid 
cycle number 9120. Open circle indicates the particular (unconverged) solution at this time. 


14 of 19 


American Institute of Aeronautics and Astronautics Paper 2006-3906 



Figure 10. SA solution map for NACA 0012 ease on 193 x 257 grid at a specific point in the boundary layer near x/c = 0.00694, at 
multigrid cycle number 11000. Open circle indicates the particular (converged) final solution. 
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Figure 11. Computed surface skin friction coefficient for the NACA 0012 airfoil on the 97 x 129 grid, SA model with fj,' t ^ = 0.009, with 
two different initial conditions, demonstrating converged solution non-uniqueness. 
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Figure 12. Computed surface skin friction coefficient for the NACA 0012 airfoil on various grids, SA model with 4,oo = 0 . 2104 . 



Figure 13. Computed surface skin friction coefficient for the NACA 0012 airfoil on various grids, SST model with fi' t oo = 0 . 009 . 


16 of 19 


American Institute of Aeronautics and Astronautics Paper 2006-3906 



multigrid cycle number 


Figure 14. Convergence of /i', k * , and uj* at specific location in the boundary layer on 193 x 257 grid, SST model with ^ = 0.009; 
phase-plane portraits will be plotted for the two times indicated in this figure. 



Figure 15. SST phase-plane portrait for NACA 0012 case on 193 x 257 grid at a specific point in the boundary layer near x/c = 0.00694, 
at multigrid cycle number 8500. Open circle indicates the particular (unconverged) solution at this time. 
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Figure 16. SST phase-plane portrait for NACA 0012 ease on 193 x 257 grid at a specific point in the boundary layer near x/c = 0.00694, 
at multigrid cycle number 11000. Open circle indicates the particular (converged) final solution. 



Figure 17. Computed surface skin friction coefficient for the NACA 0012 airfoil on various grids, SST model with /i ' t oo = 1.294. 
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Figure 18. A hypothetical phase-plane portrait for which the /c*-nullcline crosses the u;*-nullcline in the opposite direction than in Fig. 4, 
yielding a stable solution point rather than a saddle point. 
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